function q = ElnpY_CX(data, vars)

% E[ln p(Y|C, X)]

K = size(vars.Z, 2);
[N D] = size(data.Y);

DD = D*D;
Dlog2pi = D*log(2*pi);
invC = data.invC;
invC_vec = reshape(data.invC, DD, N);
logdetC = data.logdetC;
invCy = data.invCy;
yinvCy = data.yinvCy;
Z = vars.Z;

q = - Dlog2pi - logdetC;
q = q*ones(1, K);
for k = 1:K
    RRt_vec = reshape(vars.RRt(:, :, :, k), DD, N);

    q(:, k) = q(:, k) - sum(invC_vec.*RRt_vec, 1)' + sum(2*invCy.*vars.R(:, :, k), 2) - yinvCy;
    q(:, k) = q(:, k) - squeeze(sum(sum((vars.E(:, :, :, k).*invC), 1), 2));
end
q = sum(sum(Z.*q));
q = q*0.5;
